﻿Imports System.Runtime.InteropServices

Public Class GaussProjection

    Public Function GaussProjectionCal(ByVal Latitude As Double, ByVal Longitude As Double, Altitude As Double, Optional is3degrees As Boolean = True, Optional _Datum As Datum = Nothing) As GaussCoord

        'Public Latitude As Double ' 纬度 fb1,弧度
        'Public Longitude As Double ' 经度
        'Public Altitude As Double ' 高程


        Dim fx As Double ' x坐标
        Dim fy As Double ' y坐标
        Dim fB1 As Double = Latitude * Math.PI / 180 ' 纬度,弧度
        Dim fL1d As Double = Longitude ' 经度
        Dim zoneWidth As Integer = If(is3degrees = 0, 3, 6)
        Dim FL0 As Double = CalculateCentralMeridian(Longitude, zoneWidth) ' 中央经线
        Dim ZoneNo As Integer = CalculateZoneNumber(Longitude, zoneWidth) ' 带号

        If _Datum Is Nothing Then _Datum = CoordConsts.WGS84Datum

        Dim m As Double
        Dim L As Double ' 经度差，秒为单位
        L = (fL1d - FL0) * 3600
        m = Math.Cos(fB1) * L / 206264.806247096

        Dim m0 As Double
        Dim m2 As Double
        Dim m4 As Double
        Dim m6 As Double
        Dim m8 As Double
        Dim a0 As Double
        Dim a2 As Double
        Dim a4 As Double
        Dim a6 As Double
        Dim a8 As Double
        Dim ee As Double ' 第一偏心率e的平方
        ee = 2 * _Datum.F - _Datum.F * _Datum.F
        m0 = _Datum.A * (1 - ee)
        m2 = 3 * ee * m0 / 2
        m4 = 5 * ee * m2 / 4
        m6 = 7 * ee * m4 / 6
        m8 = 9 * ee * m6 / 8
        a0 = m0 + m2 / 2 + 3 * m4 / 8 + 5 * m6 / 16 + 35 * m8 / 128
        a2 = m2 / 2 + m4 / 2 + 15 * m6 / 32 + 7 * m8 / 16
        a4 = m4 / 8 + 3 * m6 / 16 + 7 * m8 / 32
        a6 = m6 / 32 + m8 / 16
        a8 = m8 / 128

        Dim zwxH As Double ' 子午线弧长
        zwxH = a0 * fB1 - a2 / 2 * Math.Sin(2 * fB1) + a4 / 4 * Math.Sin(4 * fB1) - a6 / 6 * Math.Sin(6 * fB1) + a8 / 8 * Math.Sin(8 * fB1)

        Dim n As Double ' 卯酉圈曲率半径
        n = _Datum.A / Math.Sqrt(1 - ee * Math.Sin(fB1) * Math.Sin(fB1))
        Dim e2 As Double ' e2为第二偏心率
        Dim yy As Double ' yy为η的平方
        Dim T As Double
        e2 = Math.Sqrt(ee / (1 - ee))
        yy = e2 * e2 * Math.Cos(fB1) * Math.Cos(fB1)
        T = Math.Tan(fB1)

        fx = zwxH + n * T * (m * m / 2 + (5 - T * T + 9 * yy + 4 * yy * yy) * m * m * m * m / 24 + (61 - 58 * T * T + T * T * T * T) * m * m * m * m * m * m / 720)
        fy = 500000 + n * m + n * m * m * m / 6 * (1 - T * T + yy) + n * m * m * m * m * m / 120 * (5 - 18 * T * T + T * T * T * T + 14 * yy - 58 * yy * T * T)

        Return New GaussCoord(fy, fx, FL0, ZoneNo)

    End Function

    Private Function CalculateCentralMeridian(ByVal longitude As Double, ByVal zoneWidth As Integer) As Double
        If zoneWidth = 3 Then
            Return Math.Floor((longitude + 1.5) / 3) * 3
        ElseIf zoneWidth = 6 Then
            Return Math.Floor(longitude / 6) * 6 + 3
        Else
            Throw New ArgumentException("无效的带宽。请选择3度带或6度带。")
        End If
    End Function

    Private Function CalculateZoneNumber(ByVal longitude As Double, ByVal zoneWidth As Integer) As Integer
        If zoneWidth = 3 Then
            Return CInt(Math.Floor((longitude + 1.5) / 3))
        ElseIf zoneWidth = 6 Then
            Return CInt(Math.Floor(longitude / 6)) + 1
        Else
            Throw New ArgumentException("无效的带宽。请选择3度带或6度带。")
        End If
    End Function


End Class
